En esta práctica exploraremos una aplicación en astronomía. Los datos los hemos tomado de aquí y buscaremos el mismo objetivo: predecir el diámetro de un asteroide en términos de ciertas características medidas por instrumentos satelitales y otros sensores.
La base de datos proviene del JPL (Jet Propulsion Laboratory) en Caltech y pueden buscar los datos crudos aquí. Nosotros usaremos una base de datos previamente procesada que contiene predictores para esta tarea.
Los datos tienen severos problemas de valores atípicos en la escala original de diametro (medido en metros). Esto es consecuencia de la regla de potencias descrita en la imagen de arriba. Para realizar la predicción consideraremos utilizar los datos en escala logaritmica.
Asimismo, algunos de los predictores podrían necesitar un tratamiento similar, pues algunos también incluyen valores atípicos y/o son atributos que sólo presentan valores positivos.
El objetivo es construir modelos predictivos, compararlos y escoger el mejor de ellos. Hemos desarrollado el código sobre una máquina con las siguientes características de sistema operativo, CPU (procesador y núcleos) y RAM respectivamente:
## [1] "macOS Monterey 12.2"
## [1] "Apple M1"
## [1] 8
## 8.6 GB
Por último, las predicciones se harán para un conjunto de prueba para el cual no tienen acceso a la variable objetivo. La métrica de evaluación será el error cuadrático medio en escala logarítmica, es decir,
\[ \mathscr{L}_V(y, \hat y ) = \frac1m \sum \left(\log y_i - \log \hat y_i \right)^2\,.\]
Por lo tanto, se adjuntará al proyecto un archivo *.csv con las predicciones en escala logarítmica: \(\log(\text{diametro})\).
En primer lugar, buscamos valores faltantes.
na_values <- c()
for (i in 1:ncol(train)) {
na_values[i] <- sum(is.na(train[, i]))
}
na_values## [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Al no haber datos faltantes, no necesitamos imputar por ningún método.
Como segundo paso, observamos que la mayoría de los atributos son numéricos, excepto por el atributo class.
count(train, class)## # A tibble: 9 × 2
## class n
## <chr> <int>
## 1 AMO 30
## 2 APO 44
## 3 ATE 12
## 4 CEN 3
## 5 IMB 66
## 6 MBA 12497
## 7 MCA 33
## 8 OMB 732
## 9 TJN 184
Es fácil ver que una gran mayoría de los asteroides observados son de la clase MBA tal que usar este atributo para nuestros modelos de predicción se vuelve inútil. Por lo tanto, redefinimos nuestros conjuntos de prueba y entrenamiento a escala logarítmica como sigue:
test_1 <- log(subset(test, select = -c(class)))
train_1 <- log(subset(train, select = -c(class)))El atributo condition_code se vuelve inútil después de convertirlo a la forma log(condition_code) pues toma valores desde \(-\infty\). Nótese, por el mismo efecto, su correlación con el diámetro: NaN.
train_cor <- as.data.frame(cor(train_1))
diameter_cor <- subset(train_cor, select = c(diameter))
diameter_cor %>% arrange(abs(diameter))## diameter
## w 0.0062
## om 0.0120
## ma 0.0210
## i 0.0962
## e -0.1650
## albedo -0.2504
## data_arc 0.3109
## n_obs_used 0.4493
## ad 0.4809
## moid 0.5217
## q 0.5419
## a 0.5680
## n -0.5680
## per 0.5680
## per_y 0.5680
## H -0.8430
## diameter 1.0000
## condition_code NaN
Redefinimos nuestros conjuntos como:
test_set <- subset(test_1, select = -c(condition_code))
train_set <- subset(train_1, select = -c(condition_code))Normalizamos todos los datos a través de nuestra receta para mejorar las predicciones de los modelos.
receta <- recipe(diameter ~ ., train_set) %>%
step_normalize(all_predictors()) %>%
prep()Todos los modelos siguen la siguiente estructura:
modelo_n, definimos el modelo a utilizar libre de parámetros;vc, dividimos nuestro conjunto de entrenamiento en bloques (por vfold_cv() o bootstrap());flujo_n, creamos el workflow con el modelo_n y nuestro nuestra receta;grid_n, definimos nuestro grid de búsqueda de (hiper-)parámetros (por grid_latin_hypercube(), grid_random() o de forma manual como secuencia);metricas_n, ejecutamos el proceso de ajuste de (hiper-)parámetros;mejor_n, seleccionamos las mejores métricas son select_best();modelo_final_n, redefinimos el modelo con los parámetros óptimos;flujo_final_n, creamos el workflow final con el modelo_final_n y nuestro nuestra receta;ajuste_final_n, a través de nuestro flujo_final_n ajustamos nuestro modelo a nuestro conjunto de entrenamiento;pred_n, sacamos nuestras predicciones del mismo conjunto de entrenamiento con nuestro modelo ajustado;p_n, graficamos un histograma para comparar nuestras predicciones con las etiquetas reales;q_n, hicimos un scatter plot de etiquetas reales contra nuestras predicciones;score[n], determinamos el coeficiente de determinación de cada una de las predicciones para escoger el mejor modelo a través de un score.# Modelo
modelo_1 <- linear_reg(mixture = 0, penalty = 0) %>%
set_engine("glmnet") %>%
set_mode("regression") %>%
fit(diameter ~ ., juice(receta))
coefs_1 <- tidy(modelo_1$fit) %>%
filter(term != "(Intercept)")Realizamos este plot para cerciorarnos que nuestra penalización \(\lambda\) converja.
Por lo observado en la gráfica anterior, podemos optar con seguridad por un rango de penalización entre 1e-01 y 1e+04.
modelo_1_regularizado <-
linear_reg(mixture = 0, penalty = tune()) %>%
set_engine("glmnet")
flujo_1 <- workflow() %>%
add_model(modelo_1_regularizado) %>%
add_recipe(receta)
bf_set_1 <-
parameters(penalty(range = c(-1, 4), trans = log10_trans()))
grid_1 <- grid_regular(bf_set_1, levels = 50)
vc <- vfold_cv(train_set)
tic()
metricas_1 <- tune_grid(
flujo_1,
resamples = vc,
grid = grid_1,
metrics = metric_set(rmse)
)
toc()## 4.502 sec elapsed
mejor_1 <- metricas_1 %>% select_best()Nótese la baja variabilidad en el error de validación tal que, después de pruebas con intervalos más grandes al definido en el subapartado anterior, \(\lambda \to 0\). Esto se debe a que anteriormente, hicimos una transformación logarítmica a nuestro conjunto de datos orignal.
modelo_final_1 <- finalize_model(modelo_1_regularizado, mejor_1)
flujo_final_1 <-
workflow() %>% add_model(modelo_final_1) %>% add_recipe(receta)
tic()
ajuste_final_1 <- flujo_final_1 %>% fit(train_set)
toc()## 0.186 sec elapsed
pred_1 <- predict(ajuste_final_1, new_data = train_set)
score[1] <-
R2_Score(y_pred = pred_1$.pred, y_true = train_set$diameter)La puntuación \(R^2\) obtenida con esta regresión es 0.94.
# Modelo
modelo_2 <- linear_reg(mixture = 1, penalty = 0) %>%
set_engine("glmnet") %>%
set_mode("regression") %>%
fit(diameter ~ ., juice(receta))
coefs_2 <- tidy(modelo_2$fit) %>%
filter(term != "(Intercept)")Realizamos este plot para cerciorarnos que nuestra penalización \(\lambda\) converja.
Por lo observado en la gráfica anterior, podemos optar con seguridad por un rango de penalización entre 1e-02 y 1e+01. Sin embargo, la convergencia no parece ser la deseada.
modelo_2_regularizado <-
linear_reg(mixture = 1, penalty = tune()) %>%
set_engine("glmnet")
flujo_2 <- workflow() %>%
add_model(modelo_2_regularizado) %>%
add_recipe(receta)
bf_set_2 <-
parameters(penalty(range = c(-2, 1), trans = log10_trans()))
grid_2 <- grid_regular(bf_set_2, levels = 50)
vc <- vfold_cv(train_set)
tic()
metricas_2 <- tune_grid(
flujo_2,
resamples = vc,
grid = grid_2,
metrics = metric_set(rmse)
)
toc()## 4.786 sec elapsed
mejor_2 <- metricas_2 %>% select_best()En este caso, se puede observar el mismo efecto del apartado anterior, pues nuestro error de validación \(\lambda \to 0\).
modelo_final_2 <- finalize_model(modelo_2_regularizado, mejor_2)
flujo_final_2 <-
workflow() %>% add_model(modelo_final_2) %>% add_recipe(receta)
tic()
ajuste_final_2 <- flujo_final_2 %>% fit(train_set)
toc()## 0.202 sec elapsed
pred_2 <- predict(ajuste_final_2, new_data = train_set)
score[2] <-
R2_Score(y_pred = pred_1$.pred, y_true = train_set$diameter)La puntuación \(R^2\) obtenida con esta regresión es 0.94.
# Modelo
modelo_3 <-
linear_reg(mixture = tune(), penalty = tune()) %>%
set_engine("glmnet")flujo_3 <- workflow() %>%
add_model(modelo_3) %>%
add_recipe(receta)
bf_set_3 <-
parameters(penalty(range = c(-2, 2), trans = log10_trans()), mixture(range = c(0, 1)))
grid_3 <- grid_regular(bf_set_3, levels = 50)
vc <- vfold_cv(train_set)
tic()
metricas_3 <- tune_grid(
flujo_3,
resamples = vc,
grid = grid_3,
metrics = metric_set(rmse),
control = control_grid(allow_par = TRUE) # Aquí habilitamos el cómputo en paralelo.
)
toc()## 185.223 sec elapsed
mejor_3 <- metricas_3 %>% select_best()En este caso, después de algunas pruebas, nos decidimos por difinir el intervalo de penalización de 1e-02 a 1e+02.
modelo_final_3 <- finalize_model(modelo_3, mejor_3)
flujo_final_3 <-
workflow() %>% add_model(modelo_final_3) %>% add_recipe(receta)
tic()
ajuste_final_3 <- flujo_final_3 %>% fit(train_set)
toc()## 0.189 sec elapsed
pred_3 <- predict(ajuste_final_3, new_data = train_set)
score[3] <-
R2_Score(y_pred = pred_3$.pred, y_true = train_set$diameter)La puntuación \(R^2\) obtenida con esta regresión es 0.96.
# Modelo
modelo_4 <-
decision_tree(cost_complexity = tune(), tree_depth = tune()) %>%
set_engine("rpart") %>%
set_mode("regression")grid_4 <- grid_regular(cost_complexity(),
tree_depth(),
levels = 5)
vc <- vfold_cv(train_set)
tic()
metricas_4 <- tune_grid(
modelo_4,
diameter ~ .,
resamples = vc,
grid = grid_4,
metrics = metric_set(rmse),
control = control_grid(allow_par = TRUE) # Aquí habilitamos el cómputo en paralelo.
)
toc()## 68.288 sec elapsed
mejor_4 <- metricas_4 %>% select_best()modelo_final_4 <- finalize_model(modelo_4, mejor_4)
flujo_final_4 <-
workflow() %>% add_model(modelo_final_4) %>% add_recipe(receta)
tic()
ajuste_final_4 <- flujo_final_4 %>% fit(train_set)
toc()## 0.684 sec elapsed
pred_4 <- predict(ajuste_final_4, new_data = train_set)
score[4] <-
R2_Score(y_pred = pred_4$.pred, y_true = train_set$diameter)La puntuación \(R^2\) obtenida con esta regresión es 0.98.
# Modelo
modelo_5 <-
rand_forest(mtry = tune(),
trees = tune(),
min_n = tune()) %>%
set_engine("ranger") %>%
set_mode("regression")vc <- bootstraps(train_set, times = 5)
flujo_5 <-
workflow() %>%
add_recipe(receta) %>%
add_model(modelo_5)
grid_5 <- parameters(
finalize(mtry(), select(train_set, -diameter)), # Utilizamos la función finalize()
# para que tardara menos y evitar un warning.
trees(),
min_n()) %>%
grid_random(5) # A diferencia de los *grid* de búsqueda de los demás modelos,
# en esta ocasión lo definimos de forma aleatoria.
tic()
metricas_5 <- tune_race_anova(
flujo_5,
resamples = vc,
grid = grid_5,
metrics = metric_set(rmse),
control = control_race(allow_par = TRUE) # Aquí habilitamos el cómputo en paralelo.
)
toc()## 361.568 sec elapsed
mejor_5 <- metricas_5 %>% select_best()modelo_final_5 <- finalize_model(modelo_5, mejor_5)
flujo_final_5 <-
workflow() %>% add_model(modelo_final_5) %>% add_recipe(receta)
tic()
ajuste_final_5 <- flujo_final_5 %>% fit(train_set)
toc()## 100.299 sec elapsed
pred_5 <- predict(ajuste_final_5, new_data = train_set)
score[5] <-
R2_Score(y_pred = pred_5$.pred, y_true = train_set$diameter)La puntuación \(R^2\) obtenida con esta regresión es 1.
Elegimos un modelo de máquinas de soporte vectorial con función de base radial.
# Modelo
modelo_6 <-
svm_rbf(cost = tune(),
rbf_sigma = tune()) %>%
set_engine("kernlab") %>%
set_mode("regression")vc <- bootstraps(train_set, times = 5)
flujo_6 <-
workflow() %>%
add_recipe(receta) %>%
add_model(modelo_6)
# Definimos nuestro grid de búsqueda como un hipercubo de tamaño n = 5.
grid_6 <- modelo_6 %>% parameters() %>% grid_latin_hypercube(size = 5)
tic()
metricas_6 <- tune_race_anova(
flujo_6,
resamples = vc,
grid = grid_6,
metrics = metric_set(rmse),
control = control_race(allow_par = TRUE) # Aquí habilitamos el cómputo en paralelo.
)
toc()## 159.674 sec elapsed
mejor_6 <- metricas_6 %>% select_best()modelo_final_6 <- finalize_model(modelo_6, mejor_6)
flujo_final_6 <-
workflow() %>% add_model(modelo_final_6) %>% add_recipe(receta)
tic()
ajuste_final_6 <- flujo_final_6 %>% fit(train_set)
toc()## 20.554 sec elapsed
pred_6 <- predict(ajuste_final_6, new_data = train_set)
score[6] <-
R2_Score(y_pred = pred_6$.pred, y_true = train_set$diameter)La puntuación \(R^2\) obtenida con esta regresión es 0.85.
# Modelo
modelo_7 <- boost_tree(
trees = 1000,
tree_depth = tune(),
min_n = tune(),
loss_reduction = tune(),
sample_size = tune(),
mtry = tune(),
learn_rate = tune(),
) %>%
set_engine("xgboost") %>%
set_mode("regression")vc <- vfold_cv(train_set, strata = diameter)
grid_7 <- grid_latin_hypercube( # Definimos nuestro grid de búsqueda como un
# hipercubo de tamaño 10.
tree_depth(),
min_n(),
loss_reduction(),
sample_size = sample_prop(),
finalize(mtry(), train_set),
learn_rate(),
size = 10
)
flujo_7 <- workflow() %>%
add_recipe(receta) %>%
add_model(modelo_7)
tic()
metricas_7 <- tune_race_anova(
flujo_7,
resamples = vc,
grid = grid_7,
metrics = metric_set(rmse),
control = control_race(allow_par = TRUE, parallel_over = "resamples")
) # Habilitamos el cómputo en paralelo solo en el remuestreo
# para reducir el tiempo de ejecución.
toc()## 209.883 sec elapsed
mejor_7 <- metricas_7 %>% select_best()modelo_final_7 <- finalize_model(modelo_7, mejor_7)
flujo_final_7 <-
workflow() %>% add_model(modelo_final_7) %>% add_recipe(receta)
tic()
ajuste_final_7 <- flujo_final_7 %>% fit(train_set)
toc()## 24.299 sec elapsed
pred_7 <- predict(ajuste_final_7, new_data = train_set)
score[7] <-
R2_Score(y_pred = pred_7$.pred, y_true = train_set$diameter)La puntuación \(R^2\) obtenida con esta regresión es 0.98.
## Warning: `pull_workflow_fit()` was deprecated in workflows 0.2.3.
## Please use `extract_fit_parsnip()` instead.
A pesar de que no fijamos una semilla para la aleatoriedad, el modelo de Bosques Aleatorios obtuvo siempre el mejor desempeño y las mejores predicciones. Esto nos conduce a pensar que no existe multicolinealidad entre los atributos. Por lo tanto, nuestras predicciones en escala logarítmica son las siguientes:
modelo_final <- ajuste_final_5
write_csv(as.data.frame(predict(modelo_final, new_data = test_set)$.pred),
'predicciones_04.csv')